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Abstract. We have applied a collocation approach to obtain the numerical solution 
to the stationary Schrodinger equation for systems of coupled oscillators. The 
dependence of the discretized Hamiltonian on scale and angle parameters is exploited 
to obtain optimal convergence to the exact results. A careful comparison with results 
taken from the literature is performed, showing the advantages of the present approach. 
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1. INTRODUCTION 

Coupled anharmonic oscillators have proved useful for modelling a wide variety of 
physical problems, such as, for example, the vibrations of polyatomic molecules (TJ [2J [3]. 
For this reason there has been great interest in the calculation of their eigenvalues and 
eigenfunctions [U SI El El El El [IHl HS1 HH H31 H21 HH HU] (and references therein). 

In this paper we propose an alternative variational method for the calculation 
of eigenvalues and eigenfunctions of coupled anharmonic oscillators. We develop the 
approach in Sec. in Sec. [3] we consider four coupled harmonic oscillators that are 
useful for testing our approach because the model is exactly solvable. In sections HI [5] 
and [6] we obtain the eigenvalues of two, three and four coupled anharmonic oscillators, 
respectively. In each case we compare our results with those obtained earlier by other 
authors using different approaches. Finally, in Sec. [7J we draw conclusions. 

2. THE METHOD 

In order to obtain accurate eigenvalues and eigenfunctions of the Schrodinger equation 
for coupled anharmonic oscillators we propose a collocation approach that allows the 
discretization of a D-dimensional region by means of a particular set of functions called 
Little Sine functions (LSF) pH QH UM M, [21] • They were proposed by Baye[22] who 
called them "first sine basis" . This basis set proved useful for the accurate variational 
treatment of the one-dimensional Schrodinger equation[T7], the representation of non- 
local operators on a uniform grid for the solution of the relativistic Salpeter equation 
|18j . and for the accurate treatment of the Helmholtz equation on arbitrary two- 
dimensional domains, both for the homogeneous [19] and inhomogeneous [20] case. 
Recently, we investigated the practical utility of four sets of LSF [21] with different 
boundary conditions that include the original set[T7] as a particular case. Here, we 
restrict ourselves to the LSF set with Dirichlet boundary conditions [17]. 

To make our discussion self-contained we outline the main features of our approach, 
starting with the set of functions which are used for discretization. A LSF is an 
approximate representation of the Dirac delta function in terms of the wave functions 
of a particle in a box of size 2L: 



where x±i x ) = 2WJl( x ^-kh). The index k takes all the integer values between —N/2 + 1 
and N/2 — 1, where N is an even integer. The LSF Sk is peaked at Xk = 2Lk/N = kh, 
h being the grid spacing and 2L the total extension of the interval where the function 
is defined. These LSF satisfy Sk(h, N,Xj) = 5kj and are orthogonal 




(1) 
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It follows from those properties of the LSF that we can approximate a function 
defined on x G (— L, L) as 

N/2-1 

J2 fWs k (h,N,x). (2) 

fc=-AT/2+l 

In a similar way we can also obtain a representation of the derivatives of a LSF as: 



dskih, N,x) ^ ^-^ dsk(h, N,x) 



dx ^— ' dx 

j 

S 1 ) 



Sj(h, N, x) 



= aj (h,N,x) (3) 



(i 2 s fc (/i, JV, x) d 2 s fc (fr, ^ x ) 



<ix 2 ^-^ <ix 2 



.7 



Sj(/i, iV, x) 



= 52 ci 3 s i( h > N >*) > ( 4 ) 

where the analytical expressions for the coefficients cjy have been given elsewhere |17j. 

Although Eq. (T5]) is not exact, we can make the error of that representation of the 
function f(x) as small as possible by simply increasing the value of N, as discussed 
in our earlier paper [17]. The effect of this approximation is the discretization of the 
continuous interval 2L into a set of N — 1 uniformly spaced points, x^. For example, 
the application of this approach to a one-dimensional eigenvalue problem results in the 
diagonalization of a (N — 1) x (N — 1) matrix. 

An appropriate basis set for a D-dimensional problem is given by the direct product 
of one-dimensional LSF that generates a uniform grid with spacing h (in some particular 
cases it may be more convenient to consider different spacing in different directions). A 
set of D integers (k\, A^, . . . , hp) completely specifies the location of a given point inside 
the hyper-volume. However, with the purpose of constructing Hamiltonian matrices it 
is convenient to identify one such point with just a single integer that takes all the 
values between 1 and (N — 1) D as shown in the Appendix Appendix A 



In this paper we only consider Hamiltonian operators of the form 

H = T(p 1 ,p 2 , ■ ■ .,p D ) + V{x u x 2 , ...,x D ) (5) 

where T is a polynomial of second degree and V is a polynomial function of 
the coordinates that is bounded from below. We assume that the eigenf unctions 
ijj n (xi, . . . , xd) are defined in a D-dimensional hypercube Q of side L and satisfy 
Dirichlet boundary conditons on the frontier dQ. 

By means of the discretization based on the direct product of LSF outlined above 
we obtain a Hamiltonian matrix of the form 

Hki...kn,k[...k' D — Tk 1 k' 1 $k2k' 2 ■ ■ ■ $k D k' D (6) 

+ • • • + o~kik[ ■ ■ ■ $k D - 1 k' D _ 1 Tk D k' D 

+ V(x kl ,...,x kn ) 8 klk ^...8 kDh > D (7) 
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where fc, and k[ range from —N/2 + 1 to N/2 — 1. Arguing as in Appendix Appendix A 



we express the 2D labels fc, and fc£ in terms of just two and K' D and obtain a 
(A — 1) D x (iV — 1) D matrix representation of the Hamiltonian operator. We expect that 
its eigenvalues and eigenvectors tend to the energies and wave functions, respectively, of 
the Hamiltonian operator as N increases. Notice that the potential part of this matrix 
is diagonal, while the kinetic one is sparse. These two features facilitate the numerical 
treatment of multidimensional problems as we will see in what follows. 

In principle we have two adjustable parameters: N and L, but we can bind them 
together by means of the variational method. As discussed elsewhere 1 17] it is convenient 
to set the optimal value of the scale parameter L in such a way that the trace of the 
Hamiltonian matrix is minimum. Since the Hamiltonian matrix is a relatively simple 
analytical function of L the calculation just indicated does not offer any dificulty and 
we obtain the optimal scale parameter as an analytic function of N: Lpms(N). Here, 
PMS stand for principle of minimal sensitivity [23] and the resulting Hamiltonian matrix 
depends only on N. The construction of the kinetic-energy matrix T is the time- 
consuming part of the process of building the Hamiltonian matrix H. However, in the 
problems discussed here the form of T depends only on N and D and is suitable for 
several models with different potential-energy functions V. We can thus take advantage 
of the fact that the calculation of the matrix V is faster because it is diagonal with only 
(N — 1) D elements. 

3. COUPLED HARMONIC OSCILLATORS 

In order to test the accuracy and rate of convergence of our approach we first consider 
a set of D coupled harmonic oscillators given by the Hamiltonian operator 

H = E {~\^ + \ x ^) + \ E E '•»•'•<•'•.' ( 8 ) 

where % = Vji. One can easily solve the Schrodinger equation and obtain the eigenvalues 
exactly: 



^n = ^V / l + AiU i + ^J (9) 

where rii = 0, 1, . . ., % = 1,2, ... ,D are the harmonic oscillator quantum numbers and 
Aj are the eigenvalues of the symmetrical matrix with elements Vij . 

The calculation is simple, we first construct the Hamiltonian matrix and choose 
the value of L that makes its trace a minimum [17]. Then we obtain the eigenvalues 
of the resulting matrix for increasing values of N. Table [T] shows the energies of the 
ground and first excited states of D = 4 harmonic oscillators with = (1 — <%)/3 
for increasing matrix dimension. We appreciate that the rate of convergence of present 
method is satisfactory for the treatment of coupled oscillators. With that set of potential 
parameters we obtain Ai = 1 and A 2 = A 3 = A 4 = — 1/3 so that the first excited state 
shown in Table [H is three-fold degenerate. 



1 
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N 4 x N A 




5 4 x 5 4 1.929802495 2.755212192 

7 4 x 7 4 1.931801216 2.748861410 

9 4 x 9 4 1.931851103 2.748382295 

ll 4 x ll 4 1.931851707 2.748350435 

13 4 x 13 4 1.931851659 2.748348376 

15 4 x 15 4 1.931851653 2.748348243 

exact 1.931851653 2.748348234 



15 3 1.4999999993138868 
19 3 1.4999999999986033 
29 3 1.5000000000000016 



Some time ago, Tymczak and Wangp4] calculated the eigenvalues of a three- 
dimensional harmonic oscillator (D = 3, Vij = 0) by means of different basis sets of 
wavelets. Table [2] shows that our results for the same model are considerably more 
accurate than theirs. 

4. TWO COUPLED ANHARMONIC OSCILLATORS 

In this section we apply our approach to two coupled anharmonic oscillators considered 
earlier by other authors. 

4.1. PULLEN-EDMONDS HAMILTONIAN 

Our first example of two-dimensional anharmonic oscillator is the so-called Pullen- 
Edmonds Hamiltonian [T5l [T6] 



Proceding as indicated above we obtain the eigenvalues of the Hamiltonian matrix, and 
Table [3] shows the first four of them for k — 1 . We clearly appreciate that the variational 
eigenvalues converge reasonably fast as N increases. 

Some time ago by Fessatidis and collaborators [TT] chose the Pullens-Edmons 
Hamiltonian (flUj) with k = 1 to test their proposed variational approach. The results 
obtained by those authors for the ground state, reported in their Table I, clearly converge 



Dimension 




1.4999927163457656 



H = - l -V 2 + l -(x 2 + y 2 )+Kx 2 y 2 . 



(10) 
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Table 3. Variational collocation approach to the Pullen-Edmonds hamiltonian for 
K = 1. 

iV Eq Ei Ei E% 

20 1.169791833 2.438995552 2.438995552 3.476809761 

30 1.169783302 2.438859138 2.438859138 3.475378334 

40 1.169783112 2.438854966 2.438854966 3.475320052 

50 1.169783105 2.438854795 2.438854795 3.475317137 

60 1.169783105 2.438854786 2.438854786 3.475316964 

70 1.169783105 2.438854785 2.438854785 3.475316952 



to a limit that is greater than the one obtained here by means of the LSF. The reason 
of the erroneous results of Fessatidis et al [11] lies not in the method used to solve the 
problem, but in the fact that the authors have resorted to an unsuitable spherically- 
symmetric basis set {<fj(r)}, where r = ^ x 2 + y 2 . Since the eigenfunctions of H depend 
on two variables, for example ip n (x,y) or ip n (r,<j)), x = rcos0, y = rsin0, then the 
spherical-symmetric basis generated by Fessatidis et al[H] is not complete and their 
eigenvalues do not converge to those of the Hamiltonian operator (flOl) . In fact, they 
obtained the eigenvalues of an effective central-field Hamiltonian operator in which the 
average 

i r 2w r 4 

x 2 y 2 d<p=— (11) 



substitutes the anisotropic part of the potential in Eq. ( flOi) : 

H = --V 2 + -r 2 + -r 4 (12) 
2 2 8 V 1 

We have verified that the ground-state eigenvalue of this operator for k — 1 is 
E = 1.1790711996155152844. 

4.2. QUARTIC OSCILLATORS 

Our second example of two coupled anharmonic oscillators is given by 

tt 2 2 Po, ^ 2 2 

H = 1 — miw,!, H 1 — m 2 uJnX 9 

2mi 2 1 1 2m 2 2 2 2 

+ A (c m x\ + c^x\ + c 22 x 2 ixf) , (13) 

that has been studied earlier by Hioe et al[9] and Chung and Chew [TO] for h = mi = 
m 2 = UJ\ = <jJ 2 = c 40 = c 4 = 1 and c 22 = 2. Table @] shows present results for 
the first three energy eigenvalues with different values of A and iV = 20. The reader 
may compare our results with those contained in Tables I and II of Ref.pjj], which 
also report the results of Hioe et al [H]. It is worth noticing that we can always take 
into accout the symmetry of the problem to decrease considerably the computational 
load. The eigenstates of the Hamiltonian shown above have definite parity and, for 
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Table 4. Energy eigenvalues of the system of two coupled quartic anharmonic 
oscillators for ft = mi = m<i = u>i = L02 = c±o = C04 = 1 and C22 = 2. We use a 
grid with N = 20 which yields 19 2 x 19 2 matrices. 



\ 

/A 




jiii 




yj.yjo 


1 D849Q8(sn(S 


9 9^88001 Q1 


3 4^41 (SfiDfSfi 


n 1 


1 1 ^01 881 98 


9 414^4D^(S1 


q 779QO9ROI 




1 47(^09^071 


3 931 4^3904 


^ 1 Q^^l 37Q7 


1 


1 7941 841 1 3 


3 83D3941 Q3 


fi 91 381 ^31 4 


1 n 

1U 


3 ^D1 91 D794 


7 ^97044439 


1 9 3QR81 (^9^ 


100 


0. 911899705 


15.86897394 


26.23624148 


5000 


25.27402386 


58.13369977 


96.21028659 


0.05 


1.084298606 


2.238800180 


3.454166056 


0.1 


1.150188125 


2.414340327 


3.772322591 


0.5 


1.476025046 


3.231453000 


5.195313648 


1 


1.724184069 


3.830323856 


6.213815078 


10 


3.301210571 


7.527043378 


12.39681556 


100 


6.911899338 


15.86897147 


26.23623988 


5000 


25.27402247 


58.13369048 


96.21028060 



example, we obtain the results for the even-even states shown Table H] by means of just 
10 2 x 10 2 matrices. Whenever possible the use of properly symmetrized basis sets of 
LSF is advisable in the case of large D and/or N. 

5. THREE COUPLED ANHARMONIC OSCILLATORS 

5.1. QUARTIC OSCILLATORS 

As an example of a three-dimensional anharmonic oscillator we choose the model studied 
by WitwitpH] some time ago: 

V(x, y, z) = i(> 2 + y 2 + z 2 ) + A (a xx x 4 + a yy y 4 + a zz z 4 

+ 2a xy x 2 y 2 + 2a yz y 2 z 2 + 2a xz x 2 z 2 ) . (14) 

His Table II shows the first eight eigenvalues for a xx = |, a yy = |, a zz = |, a xy = a xz = |, 
a yz = I and several values of A. Here we restrict ourselves to the most unfavourable 
case A = 10 6 . 

We have applied our collocation method to this problem following three approaches. 
The first approach is the one followed so far, which uses the minimization of the trace of 
the Hamiltonian matrix to generate the optimal scale L for a given grid. The results in 
Table [5] have been obtained in this straightforward way that does not take into account 
the anisotropy of the potential. In the second approach we take into account that 
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N 19 d x 19 



E 
Eh 
E 2 
E 3 
E 4 
E 5 
E 6 
E 7 
E 8 

Eg 



Table 5 

V(x,y,z) - 2 
A = 10 6 and a xx = h 

29 3 x 29 3 



Energy eigenvalues of the three dimensional anharmonic oscillator 
±(x 2 +y 2 + z 2 ) + \(a xx x 4 + a yy y 4 +a zz z 4 +2a X yX 2 y 2 +2ay Z y 2 z 2 + 2a xz x 2 z 2 ). 



i = i a - i 
13 w qq3 



39^ 



Ref.H 



h and a yz 



><yzi 

r 



169.2157495 
294.4522990 
315.2612020 
339.6044054 
436.2738801 
456.4724743 
487.7693071 
492.8611895 
509.1325800 
548.6572531 



169.2145773 
294.4365531 
315.2602658 
339.6041638 
436.1607904 
456.4654890 
487.7639023 
492.8571862 
509.1322591 
548.6517620 



169.2145661 
294.4363754 
315.2602614 
339.6041624 
436.1591660 
456.4654254 
487.7639454 
492.8570397 
509.1323064 
548.6516792 



169.23 
294.42 
315.28 
339.66 

456.46 

492.85 
509.14 



anisotropy by simply rescaling the y and z coordinates, y —> f3y and z — > jz, with two 
adjustable parameters (3 and 7. Thus, the resulting eigenvalue equation becomes 



+ (3 2 



of 2 



i[)(x,y/P,z/"/) 



dx 2 dy 2 dz 2 
— (E — V(x, y//3, z/-y))if>(x, y/'/3, z/ 7 ) • (15) 

As a result, now the trace of the Hamiltonian matrix explicitly depends upon /3, 7 and 
L, which are chosen to minimize it. Notice that this approach is equivalent to choosing 
LSF with different length scales on the three axes, say L x , L y and L z . Table [6] shows the 
convergence of the eigenvalues obtained in this way as the matrix dimension increases. 

The third approach consists of introducing additional adjustable parameters by 
means of a coordinate rotation of the form 



x' = x cos Q\ cos 6 2 — y sin 6\ — z cos Q\ sin 9 2 
y' = x sin Q\ cos 6 2 + y cos Q\ — z sin 9\ sin 6 2 
z' = x sin 62 + z cos 9 2 ■ 



(16) 



so that the trace of the Hamiltonian matrix now depends on L, (3, 7, 6\ and 9 2 . Only 
the potential part of the matrix will depend on the rotation angles because the kinetic 
energy is invariant under such transformation. Table [7] shows the convergence of the 
eigenvalues obtained in this way as the matrix dimension increases. 

Table [8] shows the values of the optimal parameters for the three methods outlined 
above. In Appendix Appendix B we outline some features of the variational method 



and suggest that the LSF exhibit a variational behavior that is different from that of a 
basis set of harmonic-oscillator eigenf unctions. 
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Table 6. Same as Table [5] using the second approach. 



N 


19 3 x 19 3 


29 3 x 29 3 


39 3 x 39 3 


Ref.pl] 


E 


169.2146979 


169.2145663 


169.2145660 


169.23 


Eh 


294.4375151 


294.4363709 


294.4363667 


294.42 


E 2 


315.2605725 


315.2602620 


315.2601985 


315.28 


E 3 


339.6047811 


339.6041637 


339.6041624 


339.66 


Ea 


436 1685952 


436 1591847 


436 1591447 

TitJ V.' . luJ 1 tit: ( 




E 5 


456.4659709 


456.4654324 


456.4654243 


456.46 


E 6 


487.7665445 


487.7638786 


487.7638732 




E 7 


492.8575850 


492.8576386 


492.8570724 


492.85 


Es 


509.1326192 


509.1323070 


509.1323014 


509.14 


Eg 


548.6570096 


548.6516923 


548.6516780 






Table 7. Same as Table [5] using the third approach 


N 


19 3 x 19 3 


29 3 x 29 3 


39 3 x 39 3 


Ref.[H] 


E 


169.2146303 


169.2145660 


169.2145660 


169.23 


Ex 


294.4368237 


294.4363675 


294.4363668 


294.42 


E 2 


315.2605587 


315.2602619 


315.2602616 


315.28 


E 3 


339.6043482 


339.6041626 


339.6041623 


339.66 


E, 


436.1613515 


436.1591490 


436.1591446 




E 5 


456.4664315 


456.4654261 


456.4654249 


456.46 


E 6 


487.7652490 


487.7638755 


487.7642896 




E 7 


492.8588268 


492.8571556 


492.8571528 


492.85 


E 8 


509.1332413 


509.1323079 


509.1323066 


509.14 


E 9 


548.6527375 


548.6516798 


548.6516781 





Table 8. Optimal parameters for the problems of Tables El [6] and [7] . 



TV 


LpMS 






7 




9i 


02 


20 


0.2922 














30 


0.3309 














40 


0.3623 














20 


0.2728 


0.91937 


0.1 


862^ 


17 






30 


0.3090 


0.91939 


0.1 


362£ 


53 






40 


0.3383 


0.91939 


O.i 


362^ 


11 






20 


0.2964 


1.01726 




1 




0.48115 


0.78540 


30 


0.3356 


1.01727 




1 




0.48152 


0.78540 


40 


0.3674 


1.01727 




1 




0.48164 


0.78540 
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Table 9. Energy eigenvalues of the system of three coupled sextic anharmonic 





oscillators 






N 




E l 


E 2 


5 3 x 5 3 


2.973116328 


5.292534159 


5.859553533 


9 3 x 9 3 


2.978379470 


5.296297359 


5.866068948 


17 3 x 17 3 


2.978302843 


5.295993128 


5.865822825 


19 3 x 19 3 


2.978302696 


5.295992510 


5.865822333 


21 3 x 21 3 


2.978302665 


5.295992375 


5.865822226 


29 3 x 29 3 


2.978302657 


5.295992339 


5.865822193 


Ref.[I2] 


2.978302 


5.295992 


5.865822 


Ref-HOI 


2.978305 


5.296000 


5.865828 



5.2. SEXTIC OSCILLATOR 

Our second example of three-dimensional anharmonic oscillator is given by the 
Hamiltonian operator 

H = - {pi +P 2 2 +pj) + - {x\ + x 2 2 + xf) 



+ 



x\ 



+ xl) 



+ xix 2 + X1X3 + x 2 x 3 . (17) 

studied Braun et al[12j and Chung and ( 'hew 10 . Table [9] shows our numerical results 
for several matrix dimensions and a variationally optimized grid scale L. We show the 
eigenvalues obtained by Braun et al [12] and Chung and Chew [10] for comparison. The 
last row of Table [9] shows the most accurate eigenvalues obtained by Chung and Chew 
[TO] with a matrix of size 17 3 x 17 3 . 

6. FOUR COUPLED SECTIC OSCILLATOR 



As an example of four-dimensional anharmonic oscillator we choose the model studied 
by Kaluza[13] and Chung and ChewfTU]: 



H 



- (pl+pl+pl+pl) + - (x\ + x\ + x 2 3 + xl) 



+ 



2 {x-^ ~\~ X2 ~\~ x^ ~\~ x^j ~\~ ^ (^1 ~t~ ^2 ~t~ Xq ~\~ x^ 



+ X1X2 + X1X3 + X1X4 + x 2 x 3 + X2X4 + X3X4 . (18) 

Table [10] compares our eigenvalues with those obtained by Kaluza[13] and Chung and 
Chew [TO]. Present results are more accurate than those of Chung and Chew [ID] for the 
matrix dimension 9 4 x 9 4 . 



Variational collocation for systems of coupled anharmonic oscillators 



11 



Table 10. Energy eigenvalues of the system of four coupled sextic anharmonic 





oscillators 






N 




E 1 


E 2 


3 4 x 3 4 


4.133363559 


6.144782201 


6.929503230 


5 4 x 5 4 


3.952498514 


6.276113935 


7.007385139 


9 4 x 9 4 


3.959409424 


6.281167988 


7.016036697 


ll 4 x ll 4 


3.959326310 


6.280902944 


7.015828402 


13 4 x 13 4 


3.959309441 


6.280850134 


7.015787290 


15 4 x 15 4 


3.959305195 


6.280836518 


7.015776655 


Ref.pTO] 


3.960086 


6.283305 


7.017863 


Ref.[T3| 


3.959304 







7. CONCLUSIONS 

In this paper we have extended the variational collocation approach developed earlier 
in one dimension |17j to the solution of the Schrddinger equation for systems of coupled 
oscillators or single oscillators in more than one dimensions. We have applied our method 
to several examples previously considered in the literature and obtained remarkably 
accurate results for all of them. In particular, we have shown how to improve the rate of 
convergence of our approach by means of nonlinear variational parameters. We underline 
some of the virtues of our approach: its application is straightforward and not limited 
to polynomial potentials; the construction of the Hamiltonian matrix does not involve 
the numerical calculation of integrals; one obtains its kinetic part for a given grid once 
and for all, and store it for applications to problems that differ in the potential-energy 
function; the potential part of the Hamiltonian matrix is diagonal; the Hamiltonian 
matrix is an analytic function of the variational parameters which greatly facilitates the 
numerical determination of the minimum of its trace. 

Appendix A. CONSTRUCTION OF THE HAMILTONIAN MATRIX 

The direct product of D LSF is labelled by D integers (ii, z 2 , . . . , id) where i m = 
1,2,..., M. In order to construct a M D x M D square matrix we have to encode the 
set (ii, 12, ■ ■ ■ , id) into just one integer Kd = 1,2,..., M D . One possibility is to use the 
following mapping: 

K D = M D - X {i x - 1) + M D ~ 2 (i 2 - 1) + . . . + i D (A.l) 

For the purspose of programming, it is useful to have the inverse relationship that 
produces the set (ii, i 2 , . . . ,io) for a given Kjj. Notice that we can easily extract i\ from 
Kd as follows: 

H M D ~ l + e 



(A.2) 
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where [x] stands for the integer part of x and < e < 1. Proceeding in the same way 
we have 



'2 



K D -M 



D-l. 



ll 



M D ~ 2 + e 

and similar obvious expressions for is, . . . ,io-i', finally, 



% D = K D - M 



D-l, 



^l 



M 



D-2, 



12-1] 



(A.3) 



M(i D ^ - 1) (A.4) 



The LSF label k ranges from 1 — N/2 to N/2 — 1 so that i = k + N/2 ranges from 
1 to M = N — 1 where M is odd because iV is even. 



Appendix B. VARIATIONAL METHOD 

In this section we outline some of the variational methods that one commonly applies to 
the approximate calculation of the eigenvalues and eigenfunctions of a given Hamiltonian 
operator H: 

H^ n = E n il> n , n = 0,1,... (B.l) 

One can insert a variational parameter a (or a set of such parameters) into a trial 
function ip by means of a unitary operator T(a>): (p(a) = Tip. Since T is unitary then 
A = {dT/da)T ] is antihermitian: A ] = -A. If <p is normalized then Tip will also be 
normalized and the variational method leads to the hypervirial theorem for the optimal 
variational function [25] : 

A (p(a)\H \p(a)) = (<p(a)\ [H, A] \<p( a )) = (B.2) 

If we have an orthonormal set of functions {<p n } we obtain a variational set exactly 
as before: {(j) n (a) = T<p n } and then apply the Rayleith-Ritz variational method that 
leads to 

N-l 

Y,[ H v(<*)-W6 ij ]c j = (B.3) 

The approximate energies Wk(a), k = 0, 1,...,N — 1 depend on the variational 
parameter a. Since Wk(a) > then it seems reasonable to obtain a = a°^ 1 such 
that 

^ ^0 (B.4) 

da a=a o Pt 

However, this procedure commonly requires considerable coputer time. For this reason 
many authors resort to 

dH 00 



da 



= (B.5) 

that clearly reduces computer time considerably but is mainly suitable for the ground 
state. 
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An alternative approach is to minimize the trace of the Hamiltonian matrix 

Q N-l p. N-l 

3=0 3=0 

that provides a balanced approximation to all the states considered in the Rayleigh-Ritz 
calculation. 

It is worth noticing that we can insert the adjustable parameters either into the trial 
function if or into the Hamiltonian operator according to {Tip H ftp J = (ip\ T^HT \<p). 

Consider a general dimensionless Hamiltonian operator of the form if(x, p) where 
x is a vector with components and p is a vector with the conjugate momenta pj such 
that [xj,pk] = iSjk, where i,j = 1,2, . . . ,D. In what follows we consider a variational 
basis set of the form 1 251 



f n (x) = v^et(C)0 n (Cx) (B.7) 



where C is a conveniently chosen square matrix that provides a set of variational 
parameters. Suppose that a is one of the adjustable parameters in the matrix C; then, 
it is not difficult to prove that the antihermitian operator A is given by 

j ldlndetC yr^S^ fr,-i C \ d /oca 

It is worth noticing that there is one antihermitian operator A for each variational 
parameter a although this equation does not reflect this fact explicitly. 

As a first simple example consider the scaling transformation given by = cc^y 
that leads to the scaling operators 

A = ^- + —Xi^-, i = 1, 2, . . . , D (B.9) 

If we rewrite these equations in terms of momenta in the coordinate representation 
pj = —id/dxj and rewrite the scaling transformation x'j = atjXj, p'j = ct~ x pj in 
terms of creation and anihilation operators aj = (xj + ipj)/>/2, a] = (xj — ip )/\/2, 
bj = (x'j + ipj) /y/2, bj = (x'j — ipj) /y/2 then we conclude that the Bogoliubov 
transformation bj = (hj — tjOjj / Jl — St = [a\ — tjOj^j / J~l — tj is equivalent to 
the well-known scaling discussed above. We clearly appreciate that the fashionable 
two-step approach [26 1 [TO] is simply the scaled Rayleigh-Ritz variational method that 
has proved useful since long ago in molecular calculations [27\. 

If the coordinate transformation is orthogonal C T = C _1 then the matrix 
A = C T dC j da is antisymmetric A T = —A, and 

^EEA,,^-^) (B.10) 

m n>m v " m 7 

If the problem is described by D independent coordinates Xi then we have D(D — l)/2 
components of the angular-momemtum operator like the one between parenthesis in 
Eq. flBTTOj) . 
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As an example we consider the Hamiltonian operator studied by Witwit 
H = I (pl +p 2 y +P 2 z ) + \ {x 2 + y 2 + z 2 ) 



2 xrx 2 
+ A (a xx x 4 + a yy y A + a zz z 4 

+ 2a xy x 2 y 2 + 2a xz x 2 z 2 + 2a yz y 2 z 2 ) (B.ll) 

In this case we may try three rotation parameters that will lead to the hypervirial 
relations |25j 

(<p\[H,L u ]\<p) = 0,u = x,y,z (B.12) 

If f is a product of harmonic-oscillator eigenfunctions the hypervirial theorems are 
satisfied by the unrotated functions. In other words, coordinate rotation does not appear 
to provide useful variational parameters because of the symmetry of this problem. This 
is obviously true if we use the prescriptions (1B.4j) . (IB. 51) and (IB. 61) . However, as shown 
in Sec 0, the rotation transformation is useful to improve present LSF calculations. 
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